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Abstract 



^ , We derive the mean field kinetic equation for the momentum distri- 

bution of Bogohubov excitations (bogolons) in a spatially uniform Bose- 
Einstein condensate (BEC), with a focus on the collision integrals. We 
use the method of Peletminksii and Yatsenko rather than the standard 
non-equilibrium Green's function formalism. This method produces three 
collision integrals Q^^, Q^^ and Q'^^ . Only Q^^ and Q^^ have been con- 
sidered by previous authors. The third collision integral C/^^ contains 
the effects of processes where one bogolon becomes three and vice versa. 
These processes are allowed because the total number of bogolons is not 
conserved. Since Q^^ is of the same order in the interaction strength as 
Q^^, we predict that it will significantly influence the dynamics of the bo- 
golon gas, especially the relaxation of the total number of bogolons to its 
equilibrium value. 

1 Introduction 

Kinetic equations are a primary tool in the theoretical description of non- 
uniform gases. They accurately describe the dynamics of dilute gases and are of 
theoretical interest because they provide a link between microscopic properties 
and macroscopic observables. A kinetic equation describes the evolution of the 
relevant distribution function for the system. Kinetic equations generally ac- 
count for the effects of collisions by including terms that are integrals over the 
distribution function. These terms are known as collision integrals and provide 
important information about the behavior of the gas. The structure of the colli- 
sion integral determines the conserved quantities that are used to parameterize 
the equilibrium distribution and construct hydrodynamic balance equations. 
Collision integrals also determine the rate of dissipation and relaxation and can 
be used to obtain expressions for the transport coefficients. 

Kinetic equations have a long history, beginning with the classical Boltzmann 
equation, which still provides an adequate description of dilute monoatomic 
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gases. Attempts to extend the Boltzmann equation, to treat the increasing ef- 
fects of quantum degeneracy as temperature is lowered, resulted in the Uehling- 
Uhlenbeck (UU) equation [Hll]. However, this equation fails to correctly de- 
scribe a Bose gas at all temperatures because its steady-state solution is a Bose- 
Einstein distribution in particle energies. We know that below the critical tem- 
perature, many-body effects modify the equilibrium distribution so that it de- 
pends on quasiparticle energies rather than particle energies. Many-body effects 
are accounted for by mean fields that break the U{1) gauge symmetry of the 
unperturbed Hamiltonian and necessitate a description in terms of quasiparti- 
cles. The quasiparticle description is obtained via a Bogoliubov transformation 
[3], so we refer to these quasiparticles as bogolons. 

Mean fields play a key role in the description of a condensed Bose gas, but 
determining consistent values for the mean fields is not straightforward. This 
problem has been discussed by many authors [H O [SI [Zl IE HI HH]- Some of 
the associated issues may be solved by incorporating non-uniformity [TTl [T^] 
or working in the canonical ensemble as done by V.V. and VI. V. Kocharovsky 
[131 E]. However, we note that the structure of the collision integrals only 
depends on the fact that the U{1) symmetry is broken and does not depend on 
the particular method of determining the mean fields. 

The first attempts [T5] to derive a kinetic equation for bosons that is valid at 
all temperatures were motivated by the desire to relate the equations of super- 
fluid hydrodynamics to a microscopic model of interacting bosons. These works 
did not derive the collision integrals, since they are not necessary to understand 
the gross features of superfluids. Collision integrals valid below the critical tem- 
perature were first presented by Kirkpatrick & Dorfmann [16] , and were derived 
by the method of non-equilibrium Green's functions [T^. They discussed two 
collision integrals, which we will call G^^ and G^^ that involve collisions be- 
tween the bogolons rather than particles. More recently. Griffin and coworkers 
I19[ 120] have used the non-equilibrium Greens function method to derive 
the kinetic equations of a non-uniform condensed Bose gas. Their approach is 
appropriate for describing the experimentally accessible Bose-Einstein conden- 
sates that are produced in trapping potentials. Still, the collision integrals that 
they derive are the same as those obtained by Kirkpatrick & Dorfmann. 

In this report, we shall derive the collision integrals for a condensed Bose 
gas using the alternative approach of Peletminskii & Yatsenko [21] [22] . This 
approach to kinetic theory has been used by a number of authors to describe 
relaxation processes in superfluids |23l [Ml [ISl [26] , although none of these deal 
specifically with dilute BECs. For simplicity, we shall assume that the sys- 
tem is spatially uniform. The methodical nature of this approach allows us to 
implement the derivation on a computer, which significantly reduces the time 
required to obtain the collision integrals. In addition to G^^ and G^^, we obtain 
a third collision integral (G^^) that has not been previously discussed. 

The collision integral G^^ describes a process in which one bogolon decays 
into three and vice versa. This process is important to the relaxation of the total 
bogolon number and is allowed because the number of bogolons is not conserved. 
In fact, we show that only G^'^ conserves bogolon number. All three collision 
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integrals still conserve total energy and momentum. We also show that though 
the number of bogolons relaxes, the total average number of particles remains 
constant. The nonconservation of bogolon number means that the BEC has an 
additional decay mode associated with the relaxation of bogolon number to an 
equilibrium value. Since previous works do not include G^^, only the effects of 
G^^ on the relaxation of total bogolon number have been considered. We believe 
that plays a dominant role in the relaxation of bogolon number and that it 
makes a significant contribution to other relaxation processes. 

We begin in Sec. [5] by deriving the mean- field kinetic equation using the PY 
method. In Sec. Owe introduce the mean field Hamiltonian for a condensed Bose 
gas and show how can be diagonalized using the Bogoliubov transformation. 
This diagonalized Hamiltonian governs the dynamics of a gas of bogolons. In 
Sec. |4]we give the collision integrals in the bogolon kinetic equation and in Sec. 
[5] we discuss some of their properties. We end in Sec. [5] with a summary of our 
results and concluding remarks. 

2 Particle Kinetic Equation 

The basis of the PY method of deriving kinetic equations is the Bogoliubov 
assumption [27] that for a system that is out of equilibrium, the relaxation to 
equilibrium can occur in several stages, where each successive stage has a smaller 
set of relevant parameters (expectation values and mean fields) describing the 
evolution. The density matrix is then expressed self-consistently as a function 
of the relevant expectation values and mean fields and the quantum Liouville 
equation is solved perturbatively. In this approach, setting terms of first order 
in the interaction to zero defines the mean fields, and terms of second order in 
the interaction give rise to the collision integrals. 

We consider a spatially uniform system of bosons of mass m that are confined 
to a rectangular box of volume V with periodic boundary conditions. We assume 
that the particles interact via a contact potential V{ri, rj) = gS^{ri — rj), where 
Ti is the displacement of the i**^ particle and g is the strength of the interaction. 

The Hamiltonian for this boson gas is 

H Xl^ki4iaki + ^^X!XlXl'^ki+k2,k3-k44i«k2"k3ak4, (1) 

ki ki k2 ks k4 

where eki — 2m ' ^ki creates a particle with momentum Kki and Ski destroys 
a particle with momentum Kki . The creation and annihilation operators satisfy 
the boson commutation relations [okuOk^] = <5ki,k2 where 5ki,k2 is the product 
of three Kronecker delta functions, one for each component of k. The summa- 
tions run over all single particle states for both positive and negative compo- 
nents of k. To simplify notation in subsequent sections, we will let X^ki ~^ Eii 
aj^^ — > a\, and Ski ^ '^i- We will keep this subscript convention for all quanti- 
ties which are function of the wavevector k^, such as ei — eki and Si^2 = <5ki,k2- 
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The full state of the system is described by the full density matrix p{t) which 
obeys the Liouville equation 

^h^^[H,p]. (2) 

We now implement the Bogoliubov assumption that, after a short time, the 
system evolution will relax to one governed by the behavior of the single particle 
reduced density function and the density matrix p(t) will be a functional of the 
single particle reduced density function 

r,j(i) = Tr[/}(07.,,]. (3) 

where 

Because we consider a spatially homogeneous system, we only need to consider 
diagonal elements 7^.^ of the more general operator ji j. After a sufficiently long 
time, the density operator can be written 

m = mt)), (5) 

where T{t) denotes a vector containing ri^i{t) for all values of i. The components 
Ti^i{t) are defined self-consistently such that 

T,4t) = Tr[p(r(t))7,,] = ( /"^''^^V ) (6) 

The Liouville equation ([2]) then takes the form 

'^^r- = dvzit) dt = (7) 

where r^^(t) denotes the (/x, matrix element of the 2x2 matrix Ti i{t) and 
the expression is summed over p and v. The equation for the single particle 
reduced probability density takes the form 



ih 



dt 



TT[p{Tm%^.H]). (8) 



Let us now combine Eqs. ([7|) and dS]) to obtain a self-consistent equation for 

^SSfrTr(p(rW)[7..,iI])''' = [H,mm- (9) 

Eq. ([9]) is the starting point of the derivation of the particle kinetic equation. 
In order to simplify notation for the remainder of this section, we will suppress 
the dependence of Ti^i{t) on time t. 
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The first step in solving Eq. © is the decomposition of the Hamiltonian H 
into a mean field Hamiltonian i?°(r) and an interaction term H^iT) so that H = 
H^iT) + H^iT). The form of the mean field Hamiltonian H^{T) is determined 
by the microscopic conservation laws and broken symmetries and will be written 
explicitly for BECs in the next section. It phase mixes components ^i^i of the 
single particle reduced density function such that 

[7.,.,i?"(r)] =^c,,,(r)7,-,. (10) 

J 

We can now rewrite Eq. © in terms of H^{T) and H^{T). It takes the form 

5]^Tr(p(r)[7.M,i?°(r)])'"'-[H°(r),p(r)] =^(r) (ii) 

where 

^(r) - - E ^Tr [mil.,, H\T)]y'' + [H\T),piT)]. (12) 

We can now solve Eq. (jlip in terms of a perturbation expansion in the interac- 
tion H^(T). 

Let us introduce an evolution in fictitious "time" s governed by the mean 
field Hamiltonian W^T). The s-evolution of F^.i is given by the equation 

^f^9r^^^^^^^^^^^^^^^^^^,^^^^^ ^^^^ 

Because H'^(T) commutes with itself, it is independent of s. Eq. (|lip can now 
be written 

or, more simply, 

,n^P^-[H"{T),p{T,s)]^^{T,s). (15) 
OS 

The solution p{T, s) of Eq. ([15]) is also the solution of Eq. ([TT]) for the case 
s = 0. 

Using standard methods [H], we can solve Eq. ([T5|) in terms of a perturba- 
tion expansion in the interaction H^{T). We shall assume that, after very long 
"time" s, phase mixing induced by H'^{T) and relaxation due to H^(T) cause 
the density operator to approach the limiting form 

po(r) = lini ;7"(s,0)/5(r)J7°t(s,0) = lim U°{0, s)p{T)U"'' {0, s) (16) 
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where C/°(s2,si) = e-^-^°(^(*))('*2-''i)/'*. In anticipation of the fact that the 
system wiU relax to equihbrium, we assume that po(r) can be written in the 
form 



po{T) = cxp 



(17) 



where = log (Tr [exp (— J2i -^ili)])- The quantities Xi are matrices that en- 
code the values of Ti^i. Then the solution to equation ([T5l) at "time" s = (and 
therefore also the solution to Eq. ([9])) can be written 

p{T) = po{T) + - / dsU°iO, sMT, s)U°\0, s). (18) 

The fact that Tr (^d{T)^j,j^ = implies that Tr {p{r)^jj) = Tr (/5o(r)7jj). 
This means that we can evaluate Ti^i by using the known density matrix pQ. To 
show that Tr (^^{T)jj_j^ = 0, notice that when we take the trace of the first 

term of Eq. (I12p after it has been multiplied by 7j we obtain gp''' — Sij which 
eliminates the summation. The two remaining terms in sum Eq. (jl2l) to zero. 

We now expand this solution to second order in H^{T). One can consider 
higher order terms, but the expansion to second order is sufficient for the case 
of a dilute BEG. The mean field Hamiltonian is defined so the first order con- 
tribution is zero, 

Tr(/)o(r)[7M,^'(r)]) =0. (19) 

This eliminates first order terms from the kinetic equation. To second order in 
the interaction, Eq. ([8|) takes the form 

(20) 

This is the kinetic equation that is generated by the PY method. In the next 
section, we choose a form of H^(r) that makes the first term vanish. The second 
term on the right hand side gives rise to the collision integrals. To evaluate it, we 
must now provide explicit the forms of H"{T) and H^{T) that are appropriate 
to our system. 



3 Mean Field Hamiltonian and Bogoliubov Trans- 
formation 

We can write the Hamiltonian for the BEC in a way that allows for the broken 
gauge symmetry and conserves the average number of particles N. We also want 
the form of to be such that the condition (fT9)) can be satisfied. A mean field 
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Hamiltonian H = H° + that satisfies all of these conditions can be written 
in general as 



= S + ^ (ci - /X + i/)ajai + — (at at - + a-iCii) 



and 



ijkl i 



(21) 



-S, (22) 



where the quantities u and A are mean fields given by 



= — 2^Tr(poaia^) = 



and 



^ = I^XlTr(poalaLi) = :^^Tr(poa_iai), 



(23) 



(24) 



and /U is the chemical potential. Note that though the average number of par- 
ticles is conserved by the total Hamiltonian H, it is not conserved by H°. 
The quantity S is a shift in the energy of the system and is chosen so that 
Tr ^po-ff = at zero temperature. Since S is not an operator, it has no effect 
on our derivation of the kinetic equation. 

The mean field Hamiltonian H° can be diagonalized by a Bogoliubov trans- 
formation. We define the operators 



bl = Mjd| + ViCL-i 



and 



bi = UiCli + ViCLi 



with 



Ui = 



where 



V2 



Ci+V- fl 

Ei 



+ 1, 



Vi = 



V2 



V — n 



Ei 



- 1 



Ei = ^{ei + v-iif-/^\ 



(25) 



(26) 



(27) 



The operators b\ and bj have the property that = 5ij and in terms of 

them, the mean field Hamiltonian is 



H° = Y^Ei^M + (£;. - - + M) + S. 



(28) 



Therefore, b\ and Sj are interpreted as creation and annihilation operators for 
excitations, which we will refer to as bogolons. 

It is well known that the energy spectrum of bosonic excitations must be 
gapless, that is, Ei must approach zero as i — >■ 0. This implies that n = z^— A. In 
the Green's function approach, this statement is known as the Hugenholtz-Pines 
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(HP) theorem [2S] and is expressed in terms of the self-energies as En — S12 = 
The self energies are not known exactly, but can be calculated at several different 
levels of approximation |20| . It is apparent that our mean fields will correspond 
to the self energies evaluated within one of these approximations. Setting yu = 
u — A requires us to give special treatment to the i = operators, since Ui 
and Vi are undefined (become infinite) when Ei — 0. This is to be expected 
in a BEC because the zero momentum state actually is special; it contains the 
condensate. Following Bogoliubov, we implement the special treatment of i = 
operators by replacing dj and dp with ^/Nq where Ng is the number of particles 
in the i = state. 

We can now evaluate the expectation values appearing in the self-consistency 
equations (|23)) and ((24|) as well as the collision integrals. Since the unperturbed 
Hamiltonian ff" is diagonal with respect to the operators S^S^, the density ma- 
trix pq has a Gaussian form. This allows us to use Wick theorem to evaluate 
higher-order expectation values of bogolon creation and annihilation operators. 
We also note that whenever zero appears in a particle operator momentum in- 
dex, aj and So must be replaced with ^/Nq before expectation values are taken. 

Using the Bogoliubov replacement and the Bogoliubov transformation, we 
find that the self-consistency equations become 

2g^o , 2.9^^' 2 , 2N/rti,\ , 23-^' 2 /onA 
'' = —^ + yZ^i''^+^i)iHf'^) + yZ^v, (29) 

i i 

and 

i i 

These equations are identical to those obtained for the self-energies En and E12 
evaluated in the Hartree-Fock-Bogoliubov (HFB) approximation [71 [101 122] ■ It 
is well known that the self energies of the HFB approximation do not obey the 
Hugenholtz-Pines theorem. Following the lead of previous authors for dilute 
BECs [131 [TJj , we can choose to keep the HP theorem valid by using the Popov 
approximation to the self-consistency equations. In the Popov approximation, 
A is set equal to For dilute BECs, this approximation has been shown 

to give good agreement with experiment as long as the temperature is low 
compared to the critical temperature T^, namely T < Q.QTc- Much work has 
been done on more general approximations of the self-energies [H [501 1301 [33 [32] 
and more sophisticated approximations that work well at higher temperatures 
are available. It is important to note, however, that the derivation of the collision 
integrals does not depend on the method used to obtain the mean field values. 

4 Derivation of Collision Integrals 

In this section, we turn our attention towards calculation of the collision inte- 
grals. The bogolon expectation values lj)\bi) along with A'o will form a closed 
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set of kinetic variables. In order to evaluate the kinetic equations for these vari- 
ables, we first note that the time derivative of the bogolon expectation values 
can be written 

The absence of any time derivatives on m and Vi is because their time derivatives 

can be written as = i;^^ and ^ = where 0i = ^ cosh'^ (^"^T^) ■ 

The coefficient of ^ in Eq. |3I]can be expressed in terms of bogolon occupation 
numbers and is found to vanish identically. 
Let us introduce the notation 

M = (blk) = (h-h\) = 1 + M, (32) 

N^ = (a\a^) F, = (a.aj) = 1 + iV, (33) 

and 

^^^\ ((«!«-,) + {d-^d^)) ■ (34) 

The calculation of is broken up into the calculation of and This 
provides us with more information concerning the evolution of the particle oc- 
cupation numbers Ni than calculating directly. 

We can now use Eq. and the bogolon operators to evaluate the right 
hand side of Eq. (|20|) . To accomplish this, all particle operators except a J and 
0,0 are expressed in terms of bogolon operators and the s-evolution of the bo- 
golon operators is resolved. Commutators are then evaluated and the remaining 
operators are put in normal order. Finally, the expectation values are evaluated 
using Wick expansions. This process generates a very large number of individual 
terms. The labor of condensing these terms to a manageable size was greatly 
reduced through the use of a custom computer algebra code that we developed. 

The kinetic equation for the expectation value Ni (the particle distribution) 
can be condensed to the form 

^ = CPW} + Cf{Af} + Cf'Wl (35) 

where CP, Cf^ and Cf^ are collision integrals given explicitly in appendix lAl The 
kinetic equation for the value = ^ {dld^-i + o,-idij has a similar structure 
but with different collision integrals. It can be written as 

^ = + Vf{N} + T^{M}, (36) 

where expressions for T)]^ , and can also be found in appendix Rl 

To obtain an equation for the evolution of A/i (the bogolon distribution) , we 
use Eq. (PT|) . (|35l) and ([55)1 . Putting these together results in a closed kinetic 
equation for the bogolon distribution A/i , 

^ = g^iU} + grw} + gf'W} (37) 
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where the bogolon colhsion integrals are given by 



2.3 



2 



2.3 

(38) 

Qf^W} = T^Y1 '^.+2,3+4'5(-B,+^2--B3-£;4)(W,^|3,4)'(-^.-^2AA3AA4-MAA2 J-4) 
2,3.4 

(39) 

and 



ep{AA} - ^^Y^ 5^,2+^+i5{E^ -E2^E3- Ei){Wll^^^^f{F,M2M3Mi- M^F2:F3Fi) 



2,3,4 

2 



' nV^'^ S,+2+3ASiE, + E2 + E3- E^){Wll^^^,)\F,F2T3Ni-KN2M3Ti). 
2.3,4 

(40) 

The weighting functions are given in terms of Ui and Vi by 

^1,2,3 = '"l"2U3 - U1V2U3 - V1U2U3 + U1V2V3 + V1U2V3 - Wlf2^'3, (41) 
^1^,2,3,4 = U1U2U3U4 + UiV2U3Vi + U1W2W3M4 + WlM2U3l'4 + V1U2V3U4 + ViV2V3Vi 



(42) 



and 



^1^2 34 = "1^21*3^^4 + U1U2W3U4 + U1V2U3U4 + ViV2V3Ui + V1V2U3V4 + V1U2V3V4. 

(43) 

Each of these weighting functions has specific symmetry with respect to in- 
terchanges of its indices that is shared by its collision integral. The collision 
integrals Q^^ and G^^ are identical to those generated by the Green's function 
approachpjj. The collision integral G^^ appears as a natural result of our cal- 
culation using the PY approach. 

If one sets A = 0, the weighting functions W^'j ^ and W'^"^ ^ ; equal 1 while the 
weighting function j f. ; vanishes. The self-consistency equations with A = 
imply that A^o also equals zero. What remains of the kinetic equation when 
A = is just the Uehling-Uhlenbeck equation. Therefore, the kinetic equation 
appropriate for an ideal Bose gas is recovered above the critical temperature. 

The three basic processes available to bogolons are represented by the three 
collision integrals. In G^^, one bogolon decays into two bogolons. In G^^, two 
bogolons collide elastically with each other. The collision integral G^^ describes 
one bogolon decaying into three and its inverse process. 
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5 Properties of the ColUsion Integrals 



The kinetic equation conserves the total bogolon energy given by ^/-EiA/i and 
total bogolon momentum given by '^/iJ^i- In fact, each collision integral con- 
serves energy and momentum separately. This can by shown by evaluating the 
sums 

Y^mFW} = E'^^^Pi-^} = E'^^^'HA^} = (44) 



and 



This calculation relies on the symmetry properties of the collision integral and 
the weighting functions to rename summation indices. 

One can show that the total bogolon number given by "^i'-^i ^'^^ 
served. This follows from the fact that 

J2'gl'W}^0 and E'^''{-^}^0- (46) 

i i 

However, it can still be shown that 

E'^P{-^} = (47) 

i 

meaning that Q^^ acting alone conserves bogolon number, as we would expect 
from a two-body elastic collision. 

In each of the three collision processes, the rate of the process depends upon 
the bogolon momentum distribution A/i through the factors such as J-iJ-2.N'zM4^— 
NiN2^-i^i- The evolution of the distribution function proceeds until the rate 
of each process is balanced by the rate of its inverse process and equilibrium 
is achieved. With generality, we only consider equilibrium states of zero total 
momentum. A quick analysis of the collision integrals shows that equilibrium 
distribution must be a Bose-Einstein distribution in the bogolon energies. 

< = ^^uw^ry—i (48) 

The bogolons have zero chemical potential since the total number of bogolons 
is not conserved. The effects of the particle's chemical potential are included 
through the bogolon energy spectrum Ei given by Eq. (1271) . It is easy to verify 
that is a steady state solution of the kinetic equation by substitution into 
Eq. dSIl). 

It is interesting to use the kinetic equation to evaluate the time rate of 
change of the total particle number given by ^ = X^i (^l^-j)- This sum can 
be evaluate in two parts as 

^-^+E^(«1-^)- (49) 
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Using Eq. ([55]) . we can show that 



1,2,3 

(50) 

where 

W^l.2,3 = U1U2U3 + U1U2V3 + U1V2U3 + V1U2V3 + V1V2U3 + V1V2V3 (51) 

As expected, this is non-zero because the Q^^ colhsion integral contains colhsions 
with the condensate. From the fact that the total average number of particles is 
conserved, we should expect that = —J^i'lii^l'^i)- However, to show this 
from Eq. ([55)1 . we must adopt a convention for the Bogoliubov transformation 
at zero momentum. In all we have done so far, the Bogoliubov transformation 
does not act on the zero momentum operators and the particle operators aj 
and clq are treated separately. This is equivalent to defining uq = 1 and vq = 
and taking Ao = J^a = Nq. Using these substitutions to compute from Eq. 
we obtain 

^ = E''^1.2+3^(£;i - i?2 - E3)Wil,W,%{M2Af3-AflT2T3). 

1,2,3 

(52) 

This explicitly shows that = — 'Sji^l^i) ^^'^ therefore the total average 
particle number is conserved. 



6 Conclusions 

We have derived a kinetic equation for the momentum distribution of Bogoli- 
ubov excitations (bogolons) in a spatially uniform condensed Bose gas using 
the approach of Pcletminskii and Yatsenko. The kinetic equation contains three 
collision integrals, G^^ , G^^ , and G^^ that describe the bogon processes that 
lead to decay to equilibrium of the BEC. The kinetic equation conserves total 
particle number, total energy and total momentum, but does not conserve total 
bogolon number. 

The collision integral G^^ has not been analyzed by previous authors. We 
can show (to appear elsewhere) that it has a significant influence on dissipation 
and relaxation processes in the gas. Since G^^ describes a process that does not 
conserve bogolon number, we find that its most significant influence is on the 
relaxation of the total bogolon number to its equilibrium value. However, this 
collision integral also has a measurable influence on other dissipative properties 
of the Bose-condensed gas, such as relaxation rates and transport coefficients. 
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A Derivation of Collision Operators 

In this appendix, we give explicit forms for the collision integrals C and T> found 
in Eqs. ([55)) and ([55)1 . To generate the kinetic equations for the quantities Ni 
and Aj, we use Eq. pO| to get 



and 



~dr-^^ 1^ 1^ ^1+2,3+4^5+6,7+8 
1,2,3,4 5,6,7,8 

(53) 

X / (is([a|a^a3a4, [aj(s)ai(s),4(s)ae(s)a7(s)a8(s)]) 



-dT -8h^ ^ ^ <5i+2,3+4<55+6,7+8 
1,2,3,4 5,6,7,8 



X / ds([a|4a3a4, [al(s)a^j(s),a5(s)4(s)a7(s)a8(s)]) 



(54) 

Jfc2T72 Yl '^1+2,3+4(55+6,7+8 
1,2,3,4 5,6,7,8 



X / ds{[a\dld3a4,[a-i{s)di{s),dl{s)al{s)a7{s)as{s)]) 



The calculation of the commutations and traces in Eqs. (1531) and (1541) is quite 
arduous. To speed up the process, we developed a custom symbolic algebra code 
which can store and manipulate large expressions containing both operators and 
c-number variables. The code is capable of expanding commutators, evaluat- 
ing expectation values according to Wick factorization rules, performing sums 
over Kronecker delta functions and permuting summation indices to recombine 
terms in a compact form. An output file is generated containing an algebraic 
expression that is anywhere from 70 to 500 terms long. The final simplification 
and factorization is done by hand. Following this process, we obtain the Eq. 
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dSlD with 



and 



where 



12 
2,3 



(55)' 



Cf{Af} 



+2,3+4 



S{Ei +E2-E3^ Ei)Wt 



22 



2,3,4 



2,3,4 



[t^^2,3,4(-^1-^2AA3AA4 - AA1AA2 J-3 J-4) 
- tf 2 3.4(-^-l-^-2AA_3AA_4 - AA_iAA_2-F_3^-4) 



2,3,4 



1,2+3+. 



4(5(i;i - S2 - S3 - S4)W^f ^2 3 



< [Tf2,3,4(-^lA/'2A/'3A/'4-AAiJ-2J-3J-4) 
- Tf^2^3_4(J-_iAA_2AA_3AA_4 - AA_i^_2-F_3^-4) 



J2'Sl+2+3,4S{Ei +E2 + E3- S4)W|1^2,1 



2,3,4 



[Tf"f:2,3,4(-^1-^2^3A4 - MAA2AA3 J-4) 
- tf 2 3 4(^_i J-_2^-3A/'-4 - AA_iAA_2AA_3 J--4) 



1,2,3 
B 



U1U2U3 - U1V2U3 - U1U2V3 
U1U2U3 + U1V2V3 - U1V2U3 
U1U2U3U4 + UlV2V3Ui + U1V2U3V4 
Tf^2,3,4 = UlV2U3'"4 + MlU2W3"4 + U1U2U3V4 
U1U2V3U4 + M1U2W3W4 + ''il''.'2''^3'f4 



1,2,3,4 



1,2,3,4 



(56) 



(57) 



(58) 
(59) 
(60) 
(61) 
(62) 



and T is T with each u and v interchanged. 
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For Ai we obtain Eq. ([5^ with 



2,3 



(63) 



2,3 

X {FiN2Nz -MJ-aJ-g +^-iAA_2AA_3 - AA_iJ-_2-^-3), 

^?'{-^} =7^tf + i?2 - i?3 - i?4)<5l+2,3+4W^f2,3,4f^^2,3,4 

2,3,4 

(64) 

and 

Vf{N} = ~ ^Y.'S{E, -E2-E3- E,)Si,2+3+4W^i,^,n?,2.,3A 
2,3,4 

X {FxN2N-iNi - M J-2 ^-3 + F-xN-2N--iN-i - AA_i J-_2 J"_3^_4) 

+ + i?2 + i?3 - i?4)<5l+2+3,4W^|,3,2,lf^f 2,3,4 

2,3,4 

X {TiT2J'3Mi - 7Vi7V2A/'3 -F4 + ^-1^-2^-3^^-^ - A/'_i7VL2A/'_3-7^-4) 

(65) 

where 

ri;i^2,3 = U1W2W3 + W1W2W3 - W1W2W3 + Wl"2W3 + M1U2W3 - '"1M2W3 (66) 
f^f 2,3 = V1V2U3 - U1V2V3 + U1U2V3 + U1V2U3 + W1M2W3 - V1U2U3 (67) 
^M,3,4 U1W2W3W4 + MlM2W3"4 + U1M2U3U4 ~ W1W2W3U4 — WlM2U3'"4 — ViV2V^l^) 

^1^2,3.4 ~ ViV2U3Ui + UiM2W3"4 + V1U2U3V4 — U1U2V3V4 — U1V2U3V4 ~ UiV2V^Qii) 
2,3.4 = M1W2M3W4 + W1W2W3W4 + U1W2W3W4 - W1U2W3U4 - W1U2W3U4 - WiU2W^IQ) 

When these two resuhs are combined using Eq. (jST]), we obtain Eq. ([37|) . 
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